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I. ABSTRACT 

We report the phase coexistence properties of polar- 
izable Stockmayer fluids of reduced permanent dipolcs 
|m.Q | = 1.0 and 2.0 and reduced polarizabilities a* = 
0.00, 0.03, and 0.06, calculated by a series of grand 
canonical Monte Carlo simulations with the histogram 
reweighting method. In the histogram reweighting 
method, the distributions of density and energy calcu- 
lated in Grand Canonical Monte Carlo simulations are 
stored in histograms and analyzed to construct the grand 
canonical partition function of the system. All thermo- 
dynamic properties are calculated from the grand parti- 
tion function. The results are compared with Wertheim's 
renormalization perturbation theory. Deviations be- 
tween theory and simulation results for the coexistence 
envelope are near 2% for the lower dipole moment and 
10 % for the higher dipole moment we studied. 



II. INTRODUCTION 

Stockmayer fluidsE] have long been studied as models 
for fluids with permanent dipoles, such as water, ammo- 
nia, or methyl chloride. Thermodynamic properties for 
these fluids have been calculated by theory and simu- 
lations. Attempts have been made to model real dipo- 
lar fluids by fitting potential parameters of Stockmayer 
fluids. Rowlinsono fitted experimental second virial co- 
efficients to theoretically calculated ones to find Stock- 
mayer potential parameters e, a, and \m\ for some dipo- 
lar fluids. Van Leeuwentl fitted experimental coexistence 
curves to results from computer simulations. Agreement 
between the Stockmayer potential parameters calculated 
from second virial coefficients and phase coexistence data 
is only qualitative. One of the reasons why the agreement 
is not quantitative is that fitting of the second virial co- 
efficient gives the parameters for the fluid at the limit of 
zero density while fitting of the coexistence curve gives 
the parameters for the dense fluid. The interaction of 
dipolar molecules can significantly change depending on 
the density or temperature due to the redistribution of 
electron density within a molecule in response to changes 
in the molecular environment (electrostatic induction ef- 
fect). It is essential to account for this effect when phase 
coexistence properties of dipolar fluids are calculated be- 
cause electrostatic induction is much stronger in the liq- 
uid phase than in the gas phase, and molecular interac- 
tions cannot be accurately modeled by the same (non- 



polarizable) model and parameters for both phases. One 
way to include the electrostatic induction effect on a 
model of polar fluids is to introduce polarizability. 

WertheimQ has studied the effect of polarizability on 
thermodynamic properties by using a graph-theoretical 
approach. His renormalization perturbation theoryQ was 
then extended tO|-jmxtures by Venkatasubramanian et 
al.cl. Patey et alJD performed Monte Carlo simulations 
to test the free energy calculation by Wertheim's theory 
for hard spheres with moderately large reduced dipoles 
| mg | = |mo|/v kTd 3 =1.0 (where tuq is the perma- 
nent dipole moment) and reduced mean polarizabilities 
of a* = a/d 3 = 0.03,0.06,0.1 (where k is Boltzmann 
constant, T is the absolute temperature, d is the hard 
sphere diameter and a is polarizability). Venkatasubra- 
manian et al. compared theoretically calculated coexis- 
tence properties with experiments for Stockmayer fluids 
with reduced dipoles of |m.Q | = |mo|/V e<r 3 ~ 1.0 and 
reduced polarizabilities of a* = a/a 3 ~ 0.06 (e and 
a are the Lennard- Jones parameters). The results of 
these studies agree reasonably well with Wertheim's the- 
ory. Veselya performed molecular dynamics simulations 
of polarizable Stockmayer fluids and calculated the effect 
of polarizability on thermodynamic properties such as in- 
ternaLenergy or pressure. Smit et al.Q and van Leeuwen 
et al.EHl performed Gibbs ensemble Monte Carlo simula- 
tions of coexistence properties for non-polarizable Stock- 
mayer fluids. However, simulation studies of the coex- 
istence properties of polarizable Stockmayer fluids have 
not been previously published to our knowledge. 

In this paper, we present results from Grand Canoni- 
cal Monte Carlo (GCMC) simulations of Stockmayer flu- 
ids with and without polarizability for the vapor- liq- 
uid phase coexistence properties. The results are com- 
pared to Jnhc renormalization perturbation theory by 
WcrthcimDEl. We study the models with reduced dipolcs 
of |mo| = 1.0 and 2.0 and reduced polarizabilities of 
a* =0.00, 0.03, and 0.06. Examples of estimates of the 
parameters for real fluids are |mj$| = 1-84 and a* =0.08 
for watftrtl and |mo| = 1-03 and a* =0.06 for methyl 
chlorideu. Since applications of the histogram reweight- 
ing method to phase coexistence of fluids have only re- 
cently started appearing^, we begin this paper by dis- 
cussing the principle of the method and issues related 
to its practical application to predict phase coexistence 
curves at temperatures significantly below the critical 
point. 

Deviations between theory and simulation results for 
the coexistence envelope are near 2% for the lower dipole 
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moment and 10 % for the higher dipole moment studied. 

III. GCMC HISTOGRAM RE WEIGHTING 
METHOD 

Determination of phase coexistence by Monte Carlo 
simulation requires either implicit or explicit calculation 
of the free energy of a system. The Gibbs ensemble 
methoclO is used for determination of phase equilibrium 
by implicitly minimizing the total free energy of the sys- 
tem, which is separated into two phases. Non-Boltzrasum 
sampling methods (such as thermodynamic scaling)II3't2l 
and the test particle insertion methodo are among the 
methods for explicit calculatioijuof free energy of the sys- 
tem. Ferrenberg and Swendser£3 proposed the use of the 
distribution of energy and density calculated in grand 
canonical Monte Carlo (GCMC) simulations. We refer 
to this method as "the histogram reweighting method." 
The method has rarely been applied to continuous-space 
fluids, with the exception qf-a recent study by Wilding 
for the Lennard- Jones fluidlia. We chose the histogram 
reweighting method for this study because we found it to 
be computationally more efficient than the other avail- 
able methods for our systems. 

One of the attractive features of the histogram 
reweighting method is that it can be used to construct 
the grand canonical partition function, which in turn can 
be used to obtain all thermodynamic properties, includ- 
ing the free energy or coexistence properties. Moreover, 
a single simulation at a given chemical potential \x and 
temperature T can give the thermodynamic properties at 
a range of p! and T' , by virtue of the scaling properties 
in the variables (see section on "Theoretical Basis" be- 
low). In calculating coexistence properties, it is not nec- 
essary to observe phase coexistence in a single GCMC 
simulation, since they are calculated by analyzing the 
grand canonical partition function constructed by com- 
bining the histograms. GCMC simulation is appropriate 
for calculating thermodynamic properties for a range of 
densities for molecular fluids, because it does not require 
computationally expensive volume changes. The polar- 
izable Stockmayer potential does not scale simply with 
volume. 



A. Theoretical Basis 

The basic concept behind the histogram reweighting 
method of Ferrenberg and Swendserila for GCMC simu- 
lations is reviewed here. The grand canonical partition 
function of a system with chemical potential p, volume 
V, and temperature T is written as 

S( M , V,T) = J2Y1 exp[(N/3fj,) - 0U N ]CI(N, V, U N ) (1) 

N U N 



where Vt{N, V, Un) is the number of microstates for num- 
ber of particles N, volume V and energy Un; the nota- 
tion Un emphasizes that the energy level depends on the 
number of particles. We define the chemical potential fi 
by 

0(1 = In z (2) 

where z is the activity. j3 is the inverse temperature 
(f3 = 1/kT) , where k is Boltzmann's constant. The 

denotes the summation over N from to infinity, 

JV 

and denotes the summation over all energy levels for 

u N 

each N. We perfpem GCMC simulations by Norman and 
Filinov's methodlHI and store the number of observations 
of particular N and Un in a two dimensional histogram 
fn,v,T{N, Un), which is related to the components of the 
grand canonical partition function, with a simulation- 
specific constant C, by 

WW U N )-C = exp[{N0n) - I3U N MN, V, U N ) (3) 

The thermodynamic average of a property X is calcu- 
lated by 

]T X ( N > Un)U,v,t(N, U n ) 
<X>,. V , T = N %^ — . (4) 

N U N 

Next, let us consider the grand canonical partition 
function for a different thermodynamic state with chem- 
ical potential u' and temperature T", which is written 
as 

s(//,v,r') = 

= J2J2 exp[{Np'n') - p'U N }Q(N, V, U N ) 

N U N 

-EE exp[N(P'fx' - fo) - ((3' - (3)U N ] ■ 

N U N 

■ex P [(Nf3n) - 0U N ]Q(N, V, U N ) 
= E E ^p[N(f3'n' - fa) - ((3' - (3)U N ] ■ 

N U N 

■U,v,t(N,U n )-C (5) 
The thermodynamic average of a property X is then 

< X > jU / ,V,T'= 

E u n)WU, v ,t{N, U n ) 

= N U N 

EE^wwm 

N Un 
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where W = exp[N{j3'^ - (3u) - (j3' - 0)U N ], As can be 
seen from these equations, once we know the components 
of the grand canonical partition function fn,v,T(N, Un) 
at a thermodynamic state (/i, V, T), we can construct 
a grand canonical partition function at a different ther- 
modynamic state (u\ V, T') by reweighting each com- 
ponent with exp[N(/3'(j,' - - (/?' - 0)U N ]. Then, we 
can calculate the thermodynamic properties at the new 
state (//, V, T") by averaging a thermodynamic variable 
with the reweighted components of the grand canonical 
partition function. Therefore, once we have a two dimen- 
sional histogram f^y^{N, Un) from a GCMC simulation 
at a state (/x, V, T), any thermodynamic property at any 
state (//, V, T"). The constant C, which is unknown, but 
can be determined as described below, is unimportant in 
averaging a property X as in equation |^ because it can- 
cels out in the numerator and the denominator. 



B. Determination of phase coexistence 

In the grand canonical ensemble at a sub-critical tem- 
perature, the vapor-liquid coexistence can be determined 
by finding the chemical potential that gives the same 
pressure for both phases. We define the pressure P at 
a thermodynamic state (/i, V, T) by 



P 



fi,V,T 



lnEfa,V,T) 



(7) 



If a histogram is made for a thermodynamic state (/z, V , 
T), the pressure at (//, V, T') is calculated by 



P 



■ U,v,t(N,U n ) ■ C (8) 



V 



N C7 A 



where equations [| and [7] are used. At a subcritical tem- 
perature, the pressure of each phase is thus calculated ex- 
cept for the constant C . The chemical potential // that 
gives the same pressure for both phases is then found, 
and the phase coexistence thus determined. When it is 
necessary to calculate the absolute value of the partition 
function, the simulation specific constant C needs to be 
determined by a different means. In this study, the ab- 
solute value of the constant C is fixed by calculating the 
pressure in the original simulation using the virial theo- 
rem 



P= <N> ^ T kT- 
V 



(-) 



(9) 



and equating that with the pressure calculated by equa- 
tion § 



C. Combining Results of Several Simulations 

Calculation of the thermodynamic properties at vapor- 
liquid coexistence requires the histogram over a wide 
range of N and Un- However, it is difficult to get 
ftM,V,T(N, Un) with good statistical accuracy over a wide 
range of N and Un from, one simulation. In the method 
by Norman and FilinovO or in any Metropolis schemeta, 
the configurations relevant to the given thermodynamic 
state (n, V, T) are preferentially sampled and those rele- 
vant to other states are not sampled well. Therefore, we 
perform several simulations at different thermodynamic 
states and combine the information to construct the his- 
togram fp,y,T(N, Un), which is then given with good ac- 
curacy over a wide range of N and Un- In the initial 
attempts to sample a wide range of N and Un, we per- 
form GCMC simulations at a temperature slightly above 
the critical point with various chemical potentials. The 
details are described in section 5. 

Combining two simulation results is done by fixing the 
ratio of the simulation specific constants C for two sim- 
ulations at different thermodynamic states (denoted by 
subscripts 1 and 2) by imposing 



Y exp[N(Pu (3 1 )U N \ 

U n 

■U^tA^Un)-^ 
J2 exp[N((3u - p 2 n 2 ) -((3- f3 2 )U A 



U N 



U 2 ,v.t 2 (N,U n ) ■ C 2 (10) 



at the same N and T: if the original simulations are per- 
formed at different temperatures, one or both histograms 
need to be reweighted so that the two histograms are 
compared at the same temperature. It is necessary to 
choose N and T for which the relevant configurations 
are sampled by both simulations. For convenience, we 
choose T as the average of the temperatures for the two 
original simulations and N as the number of particles at 
which the density distributions at temperature T overlap 
most. A more sophisticated way of combining multiple 
simulation results, utilizing information for a range of N 
instead oLone N, has been proposed by Ferrenberg and 
SwendsenEj. The simpler method described above was 
found to be satisfactory for our systems. 



D. Comparison to the Gibbs Ensemble method 

Both the GCMC histogram reweighting method and 
the Gibbs ensemble methoctJ can be used to obtain the 
coexistence properties of a system such as the polariz- 
able Stockmayer fluid of the present study. At first sight, 
it would seem that the Gibbs ensemble method is sim- 
pler to implement, as it requires only a single simula- 
tion per temperature at which coexistence is to be deter- 
mined, rather than the series of simulations required by 



3 



the GCMC histogram method. To determine the com- 
plete phase diagram at a series of temperatures both 
methods require several simulations. However, we have 
found that statistical uncertainties for the coexistence 
properties from the GCMC method seem to be signifi- 
cantly smaller than for a Gibbs-ensemble determination 
of the phase behavior for comparable amounts of compu- 
tational time expenditure. This statement is supported 
by the small statistical uncertainties of the coexistence 
properties calculated in Tables III and IV, which would 
have required prohibitively long Gibbs ensemble calcula- 
tions. 



IV. MODEL 

A molecule of a Stockmayer fluid is a Lennard-Jones 
interaction site with an embedded point dipole at the 
center of the molecule. For polarizable models, polari; 
ability is introduced also at the center of the moleculeQi 
Assuming that the induced dipole moment is linearly de- 
pendent on the local electric field at the center of the 
molecule, the total dipole moment of molecule i is writ- 
ten as 



m a + a- Ei 



(11) 



where is the total dipole moment, mo^ is the perma- 
nent dipole moment, Ei is the local electric field at the 
center of the molecule, and a is the pola|rizability tensor. 
The local electric field is calculated ascJ 



m 
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(12) 



where Tij is the dipole tensor for molecules i and j. If 
the vector connecting the centers of molecules i and j 
is written as and the unit tensor as Jy, the dipole 
tensor is defined as 



where 



(13) 



(14) 



The total interaction of a Stockmayer fluid is U = Ulj - 
U dp , where 



^ = £M(^) 12 -(^) 



u dp = -\ £ m * ■ Et + \ £ Ei ■ a ■ Ei 

i i 

1 \ - 



(15) 



In these equations, Ulj and Uj v denote the Lennard- 
Jones interaction and the dipole-dipole interaction, re- 
spectively, e and a are the Lennard-Jones parameters. 
Since the total dipole mement of each molecule depends 
on the dipole moments of other molecules, the energy 
is calculated by an iterative procedure. Details of the 
iterative procedure are discussed in the next section. 



V. SIMULATION DETAILS 

Throughout this study, data are presented in the re- 
duced units, denoted by the superscript (*). The units 
of energy and length are reduced by the Lennard-Jones 
parameters e and er, respectively. Reduced temperature 
is T* — kT/e. Dipole moment (m) and polarizabil- 
ity (a) are reduced by the Lennard-Jones parameters as 
m* = m/V ea 3 and a* = a/a 3 . We study Stockmayer 
fluids with permanent dipole moments |mo*| = 1.0 and 
2.0 and isotropic polarizabilities a* = 0.00, 0.03, and 
0.06. 

In GCMC simulations, at each Monte Carlo step, a 
new microstate is generated by a displacement, rotation, 
and creation or destruction of a molecule. The state thus 
generated is probabilistically accepted so that the limit- 
ing distribution of generated microstates obeys the grand 
canonical ensemble. We use 10% of the Monte Carlo 
moves for displacement, 10% for rotation, 40% for cre- 
ation, and 40% for destruction. In the energy calculation 
of polarizable models, the total interagtijm of the system 
is calculated by an iterative procedureQ'tfcl described be- 
low. 

The initial values of properties for a configuration are 
indicated by (0), and the estimates for those properties at 
the k-th iteration by (k). When a molecule is cither dis- 
placed, rotated, created or destroyed, the initial estimate 
of the electric field at the center of molecule i is 



Ei(l) 



mj(0) 



(17) 



(16) 



Then, the first estimate of the total dipole mement of 
molecule i is 

mi(l) = m ,i + a-Ei(l), (18) 
and the total electrostatic interaction is 

U dp {l) = --Y j m Q ,i-E i {l). (19) 

i 

This way, the k-th estimates of the electric field, the 
dipole mement, and the dipole-dipole interaction are 

E i {k)=Y i T ii -m j {k-l) 

rrii(k) = m 0:i + a • Ei(k), 
U dp (k) = ~^mo,i • Ei(k) (20) 
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respectively. This procedure is repeated until Ud P (k) is 
converged so that \U dp (k) - U dp {k - l)\/\U dp (k - 1)| < 
0.0001 for two consecutive iterations. Typically, the num- 
ber of iterations required for this criterion at each Monte 
Carlo move turns out to be 3 to 4 for polarizability of 
a* = 0.03 and 3 to 6 for a* = 0.06 (see tables 1 and 2). 
The initial value of dipole moment of molecule i(rrii(0)) 
at each Monte Carlo step is chosen to be the dipole 
moment before the move, except for a newly inserted 
molecule for which the dipole moment is chosen to be the 
same as the permanent dipole moment. The dipole points 
to a random direction for a newly inserted molecule. In 
order to minimize the size effect for the dipole-dipole in- 
teraction, we use the Ewald sum with 256 vectors for 
the reciprocal space termsE3 for the model with |mo*| = 
1.0 and a* = 0.00, and 514 vectors for the other models 
that we study. Approximate overall cpu-time per Monte 
Carlo step for the polarizable models is 2 times slower for 
a* = 0.03 and 2.4 times slower for a* = 0.06 than the 
cpu-time for a* = 0.00 by our code designed for polar- 
izable models with the Ewald sum. Simulations of non- 
polarizable models by our code for polarizable models 
are, in turn, approximately 4 times slower than those by 
our code specifically designed for non-polarizable mod- 
els, because the energy calculation of polarizable models 
requires the calculation of the local electric field of each 
molecule. 

When two polarizable molecules approach unphysically 
close to each other during simulation, the electrostatic 
attraction increases faster than the repulsive part of the 
Lennard-Jones potential due to the increasing magnitude 
of the total dipoles and eventually the two molecules 
overlaps. In order to avoid this effect, we set a satu- 
ration point of the total dipole equal to twice the mag- 
nitude of the permanent dipole. This treatment seems 
to be reasonable because the average of the total dipole 
moment is far smaller than twice the magnitude of the 
permanent dipole in each of the simulations we performed 
(see tables 1 and 2). For the Lennard-Jones interaction, 
the potential is cut off at half of the simulation box-Length 
and the standard long range correction is appliedu. 

During the simulations, the number of observations of 
a given number of particles and energy is counted and 
the two dimensional histogram /^y^AT, E/jv) is made. 
The grid of the histogram for energy is chosen to be 0.01 
in the reduced unit. The two dimensional histograms of 
the number of particles and the energy are analyzed to 
obtain coexistence properties for a range of temperature 
by the method described in section 2. 

In most simulations, the volume is chosen to be 
V*=216. For simulations of low densities, the volume 
is chosen to be V*=2160 so that the simulation box can 
accommodate a large enough number of particles to mea- 
sure a density difference as small as Ap* — A[l/V*} ~ 
0.0005. In combining the histograms of simulations of 
different volumes, we rescale the histogram for V*=2160 
to that for V*=216 assuming 



Q(N,V,T) = n(kN,kV,T)*, (21) 
or, in terms of the component of the histogram, 

U.yAN, U N ) = U, k v,T(kN, U kN y (22) 

which is a reasonable assumption away from the criti- 
cal point, since the logarithms of the microcanonical and 
canonical partition functions are extensive properties of 
the system. 

To fix the simulation specific constant C, we use the 
pressure calculated by the virial theorem in one of the 
simulations in the gas phase. The derivative of energy in 
terms of volume is calculated by an approximation, 

dU U(V + SV)-U(V) 

dv~ sv ■ [26) 

The SV is chosen to be about 3% of the volume V. The 
energy U(V + 8V) is calculated every time a Monte Carlo 
move is accepted. 

For each model, we first estimate the approximate loca- 
tion of the critical-paint by looking up literature values 
for similar modelsolij and by performing several short 
test simulations. Then, we perform a series of GCMC 
simulations at a temperature slightly above the critical 
temperature for various chemical potentials to sample a 
wide range of density. The reasons why we choose a tem- 
perature slightly above the critical temperature are that 
large density fluctuations near the critical point make it 
easier to sample a wide range of density in a single sim- 
ulation and that for subcritical temperatures the system 
tends to fluctuate infrequently between the gas and liquid 
densities, making sampling of both phases difficult. 

An example is illustrated in figure 1 for the model 
with \m *\ = 1.0 and a* = 0.03. At T*=1.5, we per- 
formed GCMC simulations with various chemical poten- 
tials (/j* = —7.00 ~ —1.00) to cover the density range of 
interest. The mean densities calculated from these simu- 
lations are shown in figure 1 by circles. The coexistence 
properties for a range of temperatures (T* = 1.0 ~ 1.5 
for this example) are calculated from these simulation 
results. 

In order to make sure that the simulations are sam- 
pling the configurations relevant to the phase equilibrium 
properties, we performed a new series of simulations with 
temperatures and chemical potentials that are near the 
corresponding properties at phase coexistence for tem- 
peratures lower than the temperature for the first simula- 
tions (T*=1.0, 1.1, 1.2, and 1.3). The chemical potentials 
for these simulations are chosen to be the values at co- 
existence estimated from the first simulation results for 
each temperature. The mean densities calculated from 
the new series of GCMC simulations are shown in figure 
1 by squares. 

The histograms obtained from the new series of sim- 
ulations (below the critical point) and four of the first 
simulations (near the critical point) were analyzed to 
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calculate the phase coexistence properties. The distribu- 
tions of density from the histograms used for analysis are 
shown in figure 2 and the conditions of the simulations 
are listed in table 1. The fitting of the calculated coexis- 
tence densities to the law of rectilinear diameters and a 
scaling law, assuming that dipolar fluids obey the Ising 
exponent (/3 = 0.326), is shown in figure 1 by dashed 
line. We performed five sets of simulations and analy- 
ses. The length of each GCMC simulation was 1 million 
Monte Carlo steps. The averages and the error-bars were 
calculated from the five sets of data. For error-bars, we 
use the square root of the variance of mean. Thus cal- 
culated error-bars are not, in a strict sense, statistical 
errors because there are several sources of statistical er- 
rors that propagate to the final results: for example, the 
thermodynamic states of interest here are not sampled an 
exactly equal number of times in the simulations. How- 
ever, they are still good estimates of the reliability of the 
simulations and the data analysis. 

As we can see in figure 1, the mean densities of the 
second series of simulations (squares) agree well with 
the final results (dashed line). Since the chemical po- 
tentials in the second series of simulations are based 
on the estimates from the first simulations (circles), it 
turns out that the estimates of coexistence densities from 
the first simulations are fairly accurate, considering that 
all the simulations were performed at one temperature 
(T* = 1.5). 

VI. RESULTS 

The magnitude of the average total dipole moment and 
the number of iterations necessary for the convergence of 
the energy calculation, with the criterion given in the pre- 
vious section, are listed for polarizable models along with 
the conditions used for the original GCMC simulations 
in tables 1 and 2. The number of iterations necessary 
for energy calculation turns out to be 3 to 6, depending 
on the models and the thermodynamic states. The cal- 
culated magnitude of the induced dipole is larger for the 
models with larger permanent dipoles and larger polariz- 
abilities. For the model with |m *| = 2.0 and a* = 0.06, 
the induced dipole is as large as 30 % of the permanent 
dipole at T* = 1.0 and p* ~ 0.8. 

Examples of the probability density distributions from 
the GCMC simulations are shown in figure 2 for the 
Stockmayer fluid with \m *\ = 1.0 and a* — 0.03. By 
reweighting and combining the histograms obtained from 
the simulations, we get the density distribution for any 
temperature and chemical potential. Examples of the 
density distributions at coexistence are shown in figure 3 
for the Stockmayer fluid with \m *\ = 1.0 and a* — 0.03. 
Because of the small system size of the simulations, the 
two peaks in the density distribution at coexistence over- 
lap far below the critical temperature. We calculate the 
coexistence properties up to the temperature where the 



two peaks start to overlap. The results are shown in 
figures 4 to 9, and tables 3, 4 and 5. The critical temper- 
ature increases for the models of higher dipole moment 
and higher polarizability. The heat of vaporization is cal- 
culated from the results of internal energy, pressure, and 
density at coexistence. 

We estimate the approximate values of the infinite- 
system critical temperature and density by fitting the 
simulation results to the law of rectilinear diameters and 
a scaling law, assuming that dipolar fluids obey theJsing 
exponent (/3 = 0.326). Finite-size scaling methods^ can 
be used to locate the critical point with a much higher 
accuracy than the present study, but require a series of 
simulations for different system sizes. The critical tem- 
perature increases as the polarizability increases for both 
dipoles (|mo*| = 1.0 and 2.0) that we studied (see table 
5). The effect of polarizability on the critical density is 
not pronounced. Our results seem to indicate a slight 
increase in the critical density at higher polarizabilitics, 
but the differences are comparable to the statistical un- 
certainties of the calculations. 

The calculated coexistence density, pressure, and heat 
of vaporization are compared with the renormalization 
perturbation theory by Wertheim figures 4 to 9. For the 
theoretical calculation, we follow the prescription given 
by Venkatasubramania&,et al.B and use the equation of 
state by Johnson et al.E.3 and the coefficients of the pair 
and triplet correlatioii-iunctions calculated byrElytzani- 
Stephanopoulos et al.c3 and Gubbins and Twun3 for the 
Lennard-Jones reference fluid. The agreement between 
simulation and theory is relatively good for |mj| = 1.0, 
but poor for |mjj| = 2.0. Statistical uncertainties of the 
results are quite small, confirming the computational ad- 
vantages of the histogram reweighting method. 

VII. CONCLUSIONS 

We have calculated the phase coexistence properties 
of polarizable and non-polarizable Stockmayer fluids by 
the GCMC histogram reweighting method. In the his- 
togram reweighting method, the grand canonical parti- 
tion function of the system is constructed, from which 
any thermodynamic property can be derived by statisti- 
cal thermodynamic analysis. Since the thermodynamic 
state for the grand canonical partition function can be 
continuously changed by "reweighting" the histograms, 
thermodynamic properties at thermodynamic states that 
are different from the thermodynamic state of the original 
simulation can be calculated. The statistical uncertain- 
ties of the calculated results are small, confirming the 
computational advantages of the histogram reweighting 
method over._existing methods (such as Gibbs ensemble 
Monte CarlcO). 

Our results for the coexistence properties were com- 
pared to Wertheim's renormalization perturbation the- 
ory. Differences between theoretical and simulation rc- 
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suits are within 2 % for the smaller dipole (|m.Q | = 1.0) 
but only within 10 % for the higher dipole (Im^l = 2.0) 
that we studied. 
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TABLE I. Conditions for the GCMC simulations (temper- 
ature T* , chemical potenitial /i*, and volume V) for polar- 
izable Stockmayer fluids with |mo| = 1.0 and different polar- 
izabilities (a*). The calculated average density (p*) , average 
magnitude of the total dipole (|m*|) and the average num- 
ber of iterations (kitr) necessary for calculation of the energy 
for polarizable models at each Monte Carlo step (see text) 
are also listed. The numbers in parentheses are statistical 
uncertanties, in units of the last decimal point listed. 





T* 


A** 


V 


P* 


|m*| 


fcitr 


0.03 


1.00 


-4.60 


2160 


0.0118(000) 


1.0023(00) 


3.2 




1.10 


-4.43 


2160 


0.0231(001) 


1.0039(01) 


3.3 




1.20 


-4.30 


216 


0.0410(001) 


1.0051(00) 


3.8 




1.30 


-4.17 


216 


0.0732(007) 


1.0084(01) 


3.8 




1.50 


-4.20 


216 


0.1246(018) 


1.0118(02) 


4.0 




1.50 


-4.00 


216 


0.2248(087) 


1.0199(07) 


4.0 




1.50 


-3.90 


216 


0.4118(118) 


1.0335(08) 


4.1 




1.50 


-3.72 


216 


0.5225(068) 


1.0410(07) 


4.1 




1.30 


-4.07 


216 


0.6000(060) 


1.0507(06) 


4.0 




1.20 


-4.20 


216 


0.6679(030) 


1.0572(04) 


4.0 




1.10 


-4.33 


216 


0.7306(055) 


1.0646(08) 


4.0 




1.00 


-4.50 


216 


0.7674(024) 


1.0699(05) 


4.0 


0.06 


1.00 


-4.70 


2160 


0.0106(000) 


1.0048(01) 


3.4 




1.10 


-4.53 


2160 


0.0209(001) 


1.0086(01) 


3.6 




1.20 


-4.40 


216 


0.0370(002) 


1.0109(02) 


4.5 




1.30 


-4.17 


216 


0.0806(004) 


1.0224(02) 


4.6 




1.50 


-4.10 


216 


0.2920(268) 


1.0607(48) 


5.2 




1.50 


-4.085 


216 


0.3074(282) 


1.0638(55) 


5.2 




1.40 


-4.03 


216 


0.6005(052) 


1.1213(17) 


5.7 




1.30 


-4.07 


216 


0.6666(047) 


1.1364(15) 


5.7 




1.20 


-4.20 


216 


0.7217(043) 


1.1530(16) 


5.7 




1.10 


-4.33 


216 


0.7819(046) 


1.1714(14) 


5.7 




1.00 


-4.50 


216 


0.8057(061) 


1.1844(24) 


5.7 
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TABLE II. Conditions for the GCMC simulations with the 
calculated average density, average magnitude of the total 
dipole and the number of iteration necessary for calculation of 
energy for polarizable Stockmayer fluids with |m * | = 2.0 and 
different polarizabilities. Notation is the same as for Table 1. 



a* 


T* 


P* 


V* 


P* 


|ra* 




0.03 


1.70 


-7.90 


2160 


0.0129(001) 


2.0177(004) 


3.2 




1.80 


-7.70 


2160 


0.0211(002) 


2.0264(008) 


3.2 




1.90 


-7.52 


2160 


0.0327(004) 


2.0340(010) 


3.2 




2.00 


-7.37 


216 


0.0497(010) 


2.0368(011) 


3.7 




2.10 


-7.24 


216 


0.0813(036) 


2.0513(019) 


3.6 




2.20 


-7.16 


216 


0.1305(180) 


2.0669(047) 


3.7 




2.40 


-7.00 


216 


0.2128(130) 


2.0897(037) 


3.7 




2.40 


-6.75 


216 


0.3724(158) 


2.1275(041) 


3.7 




2.20 


-6.96 


216 


0.5143(096) 


2.1641(022) 


3.8 




2.10 


-7.14 


216 


0.5723(107) 


2.1793(025) 


3.8 




2.00 


-7.27 


216 


0.6289(090) 


2.1918(020) 


3.8 




1.90 


-7.42 


216 


0.6787(037) 


2.2046(010) 


3.8 




1.80 


-7.60 


216 


0.7147(071) 


2.2158(019) 


3.8 




1.70 


-7.80 


216 


0.7540(059) 


2.2255(012) 


3.8 


0.06 


2.00 


-8.71 


2160 


0.0176(002) 


2.0443(016) 


3.4 




2.10 


-8.51 


2160 


0.0260(001) 


2.0579(007) 


3.4 




2.20 


-8.32 


2160 


0.0376(003) 


2.0766(018) 


3.5 




2.30 


-8.14 


216 


0.0555(012) 


2.0869(030) 


4.4 




o A n 
2.40 


"7 no 
- / .98 


ZlD 


0.0900(lUy) 


2.1219(114) 


A C 

4.0 




2.50 


-7.85 


216 


0.1378(136) 


2.1719(149) 


4.6 




2.70 


-7.60 


216 


0.2940(184) 


2.2762(121) 


4.8 




2.70 


-7.40 


216 


0.4318(182) 


2.3574(124) 


5.0 




2.50 


-7.70 


216 


0.5335(080) 


2.4353(063) 


5.1 




2.40 


-7.83 


216 


0.5889(059) 


2.4725(046) 


5.1 




2.30 


-8.04 


216 


0.6564(109) 


2.5199(068) 


5.2 




2.20 


-8.22 


216 


0.6927(068) 


2.5417(044) 


5.2 




2.10 


-8.41 


216 


0.7284(109) 


2.5724(076) 


5.2 




2.00 


-8.61 


216 


0.7652(084) 


2.6034(066) 


5.2 



TABLE III. Calculated coexistence properties: chemical 
potential (/!*), pressure (p*), gas phase density (p* g ), liq- 
uid phase density (pi), internal energies per molecule of gas 
phase {u* g ) and of liquid phase (it;*), heat of vaporization per 
molecule (A/i*) at different temperatures (T*) for polarizable 
Stockmayer fluids with |m *| = 1.0 and a* = 0.00, 0.03, and 
0.06. The numbers in parentheses are statistical uncertanties, 
in units of the last decimal point listed. 



* 

a 


T* 


* 


* 

P 


* 
Pa 


» 
Pi 


* 


» 


Ah* 


0.00 


1 .00 


A Af\f\l A\ 

-4.4Uyi 4 ) 


.0179(1) 


.0190(02 ) 


.754(3) 


-.276(05) 


-6.11(3) 


6.72(3 




1.05 


-4.320(3 ) 


.0253(1) 


.0208(02 ) 


.732(2) 


-.oo8(04 ) 


-5.88(2) 


6.49(2 




1.10 


A o A n { 'i \ 
-4.z4U(o ) 


.0345 (1) 


.Uo4o(Uz ) 


.708( 1) 


A n? ( HQ A 
-.4U 1 (Uo ) 


-5.65(1) 


6.20(2 




1.15 


-4.167(2) 


n A / 1 "\ 
.0400(1 ) 


f\A K 1 /'no't 
.0401(02 ) 


.682(2) 


A on / nQ \ 
-.480(03 ) 


-5.41(2) 


5.87(2 




1.20 


-4.102(2) 


.0590(1) 


.0587(01) 


.648(2) 


-.568(02) 


-5.11(2) 


5.45(2 




1.25 


A n A r /o\ 

-4.045 (2J 


.0748(1) 


n*7f c /n i \ 
.0766(01J 


.612(2) 


-.696(01 ) 


-4.79(1) 


4.95(1 




1.30 


-3.994(2) 


.0934(1) 


.1017(03) 


.570(1) 


-.891(02) 


-4.44(1) 


4.30(1 


0.03 


1.00 


-4.565(2) 


.0138(1) 


.0162(01) 


.772(1) 


-.259(02) 


-6.41(1) 


6.99(1 




1.05 


-4.470(3) 


.0203(1) 


.0219(01) 


.751(3) 


-.304(03) 


-6.19(2) 


6.79(2 




1.10 


-4.384(3) 


.0284(1) 


.0291(01) 


.729(2) 


-.349(03) 


-5.96(2) 


6.55(2 




1.15 


-4.306(4) 


.0384(1) 


.0381(02) 


.700(1) 


-.400(03) 


-5.68(1) 


6.24(1 




1.20 


-4.237(3) 


.0503(2) 


.0497(02) 


.668(3) 


-.477(03) 


-5.39(2) 


5.85(2 




1.25 


-4.176(3) 


.0645(2) 


.0650(03) 


.636(3) 


-.596(03) 


-5.10(2) 


5.40(2 




1.30 


-4.121(3) 


.0813(2) 


.0857(04) 


.601(3) 


-.770(03) 


-4.80(2) 


4.84(3 




1.35 


-4.072(2) 


.1011(2) 


.1150(05) 


.556(2) 


-1.016(04) 


-4.42(2) 


4.10(2 


0.06 


1.00 


-4.773(3) 


.0105(1) 


.0135(00) 


.797(4) 


-.241(04) 


-6.83(3) 


7.35(3 




1.05 


-4.669(3) 


.0162(1) 


.0182(01) 


.777(2) 


-.281(05) 


-6.61(2) 


7.19(2 




1.10 


-4.574(2) 


.0233(1) 


.0241(01) 


.755(3) 


-.320(05) 


-6.37(2) 


6.98(3 




1.15 


-4.487(2) 


.0320(1) 


.0315(01) 


.729(2) 


-.358(05) 


-6.11(1) 


6.73(2 




1.20 


-4.410(2) 


.0426(2) 


.0408(01) 


.701(1) 


-.411(03) 


-5.83(1) 


6.41(1 




1.25 


-4.341(2) 


.0551(2) 


.0528(01) 


.668(2) 


-.498(01) 


-5.52(2) 


5.99(2 




1.30 


-4.279(3) 


.0699(3) 


.0688(04) 


.636(2) 


-.634(03) 


-5.23(2) 


5.50(2 




1.35 


-4.224(4) 


.0874(4) 


.0908(09) 


.602(2) 


-.828(09) 


-4.92(1) 


4.91(2 




1.40 


-4.174(4) 


.1081(6) 


.1218(23) 


.555(4) 


-1.097(22) 


-4.52(3) 


4.11(4 



TABLE IV. Calculated coexistence properties for polariz- 
able Stockmayer fluids with |rao*| = 2.0 and a* = 0.00, 0.03, 
and 0.06. Notation is the same as in table 3. 



0.00 



0.03 



0.06 



T* 



1.60 
1.65 
1.70 
1.75 
1.80 
1.85 
1.90 



1.70 
1.75 
1.80 
1.85 
1.90 
1.95 
2.00 
2.05 
2.10 



2.00 
2.05 
2.10 
2.15 
2.20 
2.25 
2.30 
2.35 
2.40 



P 



Pi 



-7.177(3 
-7.078(4 
-6.986(4 
-6.900(4 
-6.822(4 
-6.750(4 
-6.683(4 



-7.893(6 
-7.783(4 
-7.679(3 
-7.582(2 
-7.492(3 
-7.407(3 
-7.329(4 
-7.255(4 
-7.187(5 



-8.695(8 
■8.581(8 
-8.471(8 
-8.367(8 
-8.270(7 
■8.178(7 
-8.092(6 
-8.010(6 
-7.934(5 



.0224(0 
.0300(1 
.0388(1 
.0490(1 
.0607(2 
.0741(3 
.0896(3 



.0233(4 
.0301(4 
.0379(4 
.0469(4 
.0572(5 
.0690(5 
.0823(5 
.0974(6 
.1146(7 



.0357(4 
.0440(4 
.0535(4 
.0641(4 
.0761(4 
.0895(4 
.1044(4 
.1211(5 
.1397(6 



.0222(01 
.0273(02 
.0335(02 
.0411(03 
.0507(04 
.0634(05 
.0811(07 



.0183(01 
.0222(01 
.0270(01 
.0327(00 
.0396(00 
.0480(01 
.0587(02 
.0730(05 
.0927(11 



.0232(02 
.0275(02 
.0325(02 
.0385(02 
.0455(03 
.0539(04 
.0644(05 
.0784(07 
.0969(12 



.726(2 
.706(2 
.682(2 
.654(2 
.626(2 
.599(2 
.569(4 



.758(4 
.741(5 
.720(6 
.697(5 
.674(4 
.650(4 
.624(4 
.593(4 
.556(4 



.761(1 
.749(1 
.732(3 
.709(4 
.686(4 
.663(3 
.639(4 
.609(5 
.573(7 



-1.16(1 
-1.28(1 
-1.40(1 
-1.53(1 
-1.67(1 
-1.86(1 
-2.13(2 



-1.11(1 
-1.22(1 
-1.34(1 
-1.46(1 
-1.58(1 
-1.71(1 
-1.87(0 
-2.10(1 
-2.44(2 



-1TT9T2 
-1.31(2 
-1.44(2 
-1.58(2 
-1.71(2 
-1.86(2 
-2.05(2 
-2.30(2 
-2.66(2 



-10.16(2 
-9.89(3 
-9.58(3 
-9.22(3 
-8.87(2 
-8.53(2 
-8.17(4 



-11.48(6 
-11.21(6 
-10.91(7 
-10.58(6 
-10.25(5 
-9.92(4 
-9.56(4 
-9.16(4 
-8.69(4 



-13.00(3 
-12.75(4 
-12.44(5 
-12.07(4 
-11.69(4 
-11.33(4 
-10.95(5 
-10.51(7 
-9 



Ah' 
~ 9797(2 
9.67(3 
9.28(4 
8.82(3 
8.30(3 
7.72(3 
6.98(4 
11.62(9 
11.30(8 
10.91 
10.48(8 
10.02(6 
9.53(5 
8.95(5 
8.21(5 
7.27(6 
13.30(4 
12.98(5 
12.57(5 
12.07(5 
11.54(5 
11.00(5 
10.36(6 
9.55 
8.53(9 
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TABLE V. Estimates of critical temperature and density 
for polarizable Stockmayer fluids as a function of the perma- 
nent dipole (|mo*|) and polarizability (a*). 



mo* 


* 

a 


T* 

' cr 


Per 


1.00 


0.00 


1.400(3) 


0.318(3) 




0.03 


1.432(6) 


0.322(5) 




0.06 


1.478(8) 


0.328(6) 


2.00 


0.00 


2.05(1) 


0.306(08) 




0.03 


2.22(1) 


0.302(09) 




0.06 


2.53(2) 


0.312(16) 



1.6 
1.4 
1.2 
1.0 t 
0.8 



-> 1 1 1 1 1 ' r 

o o p o o oo 



El 
/ 



\ 



a. 



hi 



_i i i . i . i_ 



0.0 0.2 0.4 0.6 0.8 

* 

P 



FIG. 1. Temperature and mean density of example GCMC 
simulations. The initial simulations are performed to cover a 
wide range of density at a temperature slightly higher than 
the estimated critical point (T* = 1.5) with chemical poten- 
tials of n*=-7.00, -6.00, -5.00, -4.20, -4.00, -3.90, -3.72, -3.00, 
-2.00, and -1.00 (circles, from left to right). The second series 
of simulations are performed at temperatures and chemical 
potentials for phase coexistence (squares, see also table 1), 
as calculated from the initial simulations. Dashed line is the 
fitting of the coexistence curve to the results from the second 
series of simulations. 



60 



40 



20 




0.0 0.2 0.4 0.6 0.8 

* 

P 



FIG. 2. Density distributions from the GCMC simulations 
with |m *| = 1.0 and a* = 0.03. The chemical potential and 
the temperature for each simulation is, from left to right, (/tt* = 
-4.6, T*=1.0), (-4.43, 1.1), (-4.30,1.2), (-4.17,1.3), (-4.20,1.5), 
(-4.00,1.5), (-3.90,1.5), (-3.72,1.5), (-4.07,1.3), (-4.20,1.2), 
(-4.33,1.1), (-4.50,1.0) (see also table 1). The whole density 
range of interest is sampled by these simulations. 
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FIG. 3. Density distributions at two-phase coexistence 
for the polarizable Stockmayer fluid with |m *| = 1.0 and 
a* — 0.03, calculated by reweighting the histogram ob- 
tained from the GCMC simulations (see table 1 and figure 
1). The distributions are shown for T* = 1.00 (solid line), 
1.20 (dot-dashed line), and 1.40 (dashed line). 



FIG. 5. Coexistence densities for the polarizable Stock- 
mayer fluids with | mo* | = 2.0. The notation is the same 
as in figure 4. 




0.2 0.4 0.6 0.8 
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FIG. 4. Coexistence densities for the polarizable Stock- 
mayer fluids with |m *| = 1.0. Circles, squares, triangles 
are for a* = 0.00, 0.03, and 0.06, respectively. The error bars 
are about the same or smaller than the size of the symbols. 
Solid line, dashed line, and dot-dashed line are the results of 
the Wertheim's perturbation theory for a* — 0.00, 0.03, and 
0.06 respectively. 
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FIG. 6. Coexistence pressure for the polarizable Stock- 
mayer fluids with |m *| = 1.0. The notation is the same 
as in figure 4. 
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FIG. 7. Coexistence pressure for the polarizable Stock- FIG. 9. Heat of vaporization for the polarizable Stock- 

mayer fluids with |m *| = 2.0. The notation is the same mayer fluids with |m *| = 2.0. The notation is the same 
as in figure 4. as in figure 4. 




11 



